Specific heat of quasi-2D antiferromagnetic Heisenberg models with varying 

inter-planar couplings 
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We have used the stochastic series expansion (SSE) quantum Monte Carlo (QMC) method to 
study the three-dimensional (3D) antiferromagnetic Heisenberg model on cubic lattices with in- 
plane coupling J and varying inter-plane coupling J± < J. The specific heat curves exhibit a 3D 
ordering peak as well as a broad maximum arising from short-range 2D order. For J± <C J, there is 
a clear separation of the two peaks. In the simulations, the contributions to the total specific heat 
from the ordering across and within the layers can be separated, and this enables us to study in 
detail the 3D peak around T c (which otherwise typically is dominated by statistical noise). We find 
that the peak height decreases with decreasing J±, becoming nearly linear below J± = 0.2 J. The 
relevance of these results to the lack of observed specific heat anomaly at the ordering transition of 
some quasi-2D antiferromagnets is discussed. 

PACS numbers: PACS: 75.40.Gb, 75.40.Mg, 75.10.Jm, 75.30.Ds 



I. INTRODUCTION 

Spatially anisotropic systems and dimensional 
crossovers have been issues of theoretical and experi- 
mental interest for many decades, especially in context 
of classical critical phenomenai^ In recent years, a large 
number of quasi-low dimensional, low-spin, spatially 
anisotropic materials have been synthesized and their 
properties investigated in great detail. This has led 
to a renewed interest in these issues including the 
role of enhanced quantum fluctuations Perhaps 
the most studied of these are the cuprate family of 
materials, whose parent stoichiometric compounds are 
antiferromagnetic insulators which upon doping become 
high temperature superconductors. These are layered 
compounds, where exchange coupling between the planes 
is many orders of magnitude smaller than the exchange 
coupling in the planes,-'- However, these are by no 
means the only systems where spatial anisotropy and 
dimensional crossovers are important. The list of just 
novel transition-metal oxide materials, which despite 
their low-dimensionality often develop 3-dimcnsional 
long-range order, includes several cuprates, vanadates, 
copper-germenates, pnictide oxides, manganites, etc^ 

In these materials both spatial anisotropy and 
anisotropy in spin-space can be important in the develop- 
ment of 3D order. For example, it is quite possible that in 
some cuprate families XY anisotropy plays an important 
role in bringing about long-range order, while in others it 
is the interplanar coupling which is primarily responsible 
for the transition. Here, we will focus on layered systems 
with SU(2) symmetry in spin-space. This is believed to 
be relevant to the material La2Cu04. At the finite tem- 
perature 3D transition, one expects the universality class 
for such a system to be that of classical 3D Heisenberg 
model. However, in La2Cu04 no specific heat anomaly is 
seen at the 3D transition, 10 contrary to expectations for 



the 3D classical Heisenberg model. In this paper we use 
a quantum Monte Carlo (QMC) method to verify that 
the transition in spatially anisotropic systems remains in 
the universality class of 3D classical Heisenberg model. 
Our primary goal is to clarify how the amplitude for the 
specific-heat anomaly at the transition is diminished in 
systems with weak interplanar couplings. This would 
help us predict which of the newly synthesized systems 
should show such anomalies, given the finite experimen- 
tal resolution. 

A simple way to understand the reduction in the ampli- 
tude for the specific heat anomaly, in these systems, is to 
consider the effect of preexisting short-range order at the 
transition. In a spatially anisotropic system, short range 
order in the planes can develop at temperatures much 
above the 3D ordering temperature. And, if the system is 
highly anisotropic, substantial spin-correlations can de- 
velop in the planes before the eventual 3D transition. 
This means the effective number of degrees of freedom 
involved in the 3D order is substantially reduced. Hence, 
the specific-heat anomaly must diminish. Our goal is to 
obtain a quantitative estimate for this effect. 

The rest of the paper is organized as follows. We intro- 
duce the model and the computational techniques used 
in Sec. II. The results of the simulations and the related 
discussions are presented in Sees. Ill and IV. We conclude 
in Sec. V with a summary of the results. 



II. MODELS AND SIMULATION TECHNIQUE 

We have studied the Heisenberg antiferromagnet on 
an anisotropic cubic lattice. This model is given by the 
Hamiltonian 



H — J ^ ' Si ■ Sj + J± ^ ' Sj • Sj 
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FIG. 1: Spin stiffness vs. temperature for different systems 
with the same aspect ratio. The upper (lower) panel shows 
the stiffness perpendicular (parallel) to the planes. 



(processes) with adjacent values of T at regular intervals 
(typically after every Monte Carlo step, each time at- 
tempting several hundred swaps) according to a scheme 
that maintains detailed balance in the space of the paral- 
lel simulations. This has favorable effects on the simula- 
tion dynamics, as the temperature of the SSE configura- 
tions will fluctuate across the critical temperature. More 
importantly in the case considered here, a given configu- 
ration will contribute to measured expectation values at 
several nearby temperatures, thereby reducing the over- 
all statistical errors (at the cost of introducing correla- 
tions between the errors, which is of minor significance 
here) . Implementation of tempering schemes in the con- 
text of the SSE method have been discussed in Ref. |2J. 

The thermodynamics of the 3D Heisenberg model on 
an isotropic simple cubic lattice are fairly well understood 
from both analytic and computational studies 
Recent large scale Monte Carlo studiesi 4 *^ have re- 
sulted in an accurate estimate of the critical tempera- 
ture, T c /J sa 0.946. Several approximations also exist 
for T c of the anisotropic modeliii^iSiii For weak cou- 
pling between the planes, the interplanar couplings can 
be treated in mean-field theory and lead to the relation 
T c ~ — l/ln(a)£ We are not aware of any previous cal- 
culations of the specific heat of anisotropic systems. 



III. LOCATING THE TRANSITION 
TEMPERATURE 



where J (J±) is the strength of the intra- (inter-) planar 
coupling. The first (second) summation refers to sum- 
ming over all nearest neighbors parallel (perpendicular) 
to the XY-plane. We will study the model as a function 
of the dimensionless inter-plane coupling a — Jj_/J. 

The stochastic series expansion (SSE) methodi2iS& is 
a finite-temperature QMC technique based on impor- 
tance sampling of the diagonal matrix elements of the 
density matrix e~@ H . There are no approximations be- 
yond statistical errors. Using the "operator- loop" cluster 
update^ the autocorrelation time for the system sizes we 
consider here (up to rj 3 x 10 4 spins) is at most a few 
Monte Carlo sweeps even at the critical temperaturei^i 

On the dense temperature grids that we need in 
order to study the critical region in detail, we have 
further found that the statistics of the data obtained 
can be significantly improved by the use of a temper- 
ing schemei22i2£ A standard single-process tempering 
method, where the temperature of the simulation fluc- 
tuates on a grid of pre-selected temperatures, was pre- 
viously used in a study of the isotropic 3D Heisenberg 
models Here we use parallel tempering,2£ where several 
simulations are run simultaneously on a parallel com- 
puter, using a fixed value of a and different, but closely 
spaced, values of T at and around the critical temper- 
ature. Along with the usual Monte Carlo updates, we 
attempt to swap the temperatures of SSE configurations 



We first determine the transition temperature for the 
model as a function of a. An efficient way to do this is 
by studying the scaling properties of the spin-stiffness. 
We have evaluated the spin stiffnesses both parallel to 
and perpendicular to the planes. The stiffness can be 
defined2S*S£ as the second derivative of the free energy 
with respect to a uniform twist <fi: 



d 2 F{cj>) 



(2) 



The stiffness can also be related to the fluctuations of 
the "winding number" in the simulationsi2i2L22i22i and 
hence can be estimated directly without actually includ- 
ing a twist. Since the twist can be applied parallel to or 
perpendicular to the planes, there are two different spin 
stiffnesses, p x and p z , in the anisotropic system consid- 
ered here. 

For a system of weakly coupled Heisenberg chains, 
it has been shown that estimates for various observ- 
ables for a spatially anisotropic system can depend non- 
monotonically on the system size for square lattices^ 
One can instead use rectangular lattices to more rapidly 
obtain monotonic behavior of the numerical results for 
extrapolating to the thermodynamic limit. We expect 
similar effects in the present model at a <C 1. Hence we 
have studied tetragonal lattices with L x — L y ^ L z . Lat- 
tices with an aspect ratio R = L x / L z = 4 have been used 
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to obtain the results presented here. We have chosen six 
different values o f a, of the form of a — 2~ n , n — 1, . . . , 6. 

Following Ref. 114 , we use the finite-size and tempera- 
ture dependence of the spin stiffnesses to determine the 
critical temperature^ For a fixed aspect ratio, the stiff- 
ness at T c is predicted to scale as 



L 



2-d 



(3) 



where d is the dimensionality of the system. The above 
relation implies that for the 3D Heisenberg model, on a 
plot of L^p^ as a function of T the curves for differ- 
ent system sizes will cross each other at T c . Results 
for a = 1/4 are shown in Fig. ^ The upper (lower) 
panel shows L x p x {L z p z ) versus T for four different sys- 
tem sizes. The curves indeed intersect each other almost 
at a single point. Subleading corrections are seen in the 
fact that the crossing points move slightly as the system 
size is increased. Interestingly, the behavior is opposite 
for the two stiffness constants; in the case of p x the cross- 
ings move down in temperatures, whereas the p z cross- 
ings move up. Hence, we believe that the crossings for 
the two largest system sizes bracket the true T c and we 
view them as the upper and lower bounds. From these 
results we estimate T c = 0.6160 ± 0.0005 for a = 1/4. 

Next we study the universality class of the transition. 
To this end, we consider the static magnetic susceptibil- 
ity, defined as 



X(q) 



N ^ 



,iq-(r 



dr<5f(r)fi?(0)>, (4) 



where N — L x L y is the size of the system. At the 
critical temperature, the staggered susceptibility x(Q) 
should scaled with the system length as L x ~ v , where 
Q = (7T, 7r, 7r) is the 3D ordering wave vector. For any 
non-zero value of J±, the transition is expected to be- 
long to the classical 3D Heisenberg universality class, for 
which the critical exponents are known to a high degree 
of accuracy^ The spin-spin correlation function expo- 
nent r\ w 0.037. Figure Ufa) shows a = 1/4 results for 
ln(%(Q)/i^) versus h\(L x ) at temperatures close to T c . 
Asymptotically, we expect the data to fall on a straight 
line with slope —r\ w —0.037 at T = T c and diverge up- 
ward (downward) for T < T c (T > T c ). This is indeed 
what we observe. The curves are completely consistent 
with the known value of r\ and the estimate of T c obtained 
from Fig. ^ 

We have also tested the expected scaling for T > T c . 
In the thermodynamic limit, x(Q) should diverge as i~ 7 , 
where t = \T — T c \ and 7 = v(2 — rf). For a finite system, 
finite-size scaling predicts xl(£) = Xoo(t)f[£(t)/L], with 
the correlation length diverging as £ ~ t~ v . Hence on 
a plot of XlM* 7 versus Lt u , data for different L should 
collapse onto a single curve. As shown in Fig.^b), this 
is indeed the case with our estimated T c and the known 
3D Heisenberg exponents. 

We have here discussed the determination of T c and 
checked the consistency with the expected universality 
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FIG. 2: Finite-size scaling of the staggered susceptibility at 
a = 1/4. (a) Size dependence close to T c . At T c , the data are 
expected to fall on a straight line with slope —r\ = —0.037, 
which is indicated with the dotted line, (b) Scaling plot above 
T c , using r c /J = 0.616 and the 3D classical Heisenberg expo- 
nents ?? = 0.037 and v = 0.711. 



class for a — 1/4. Using the spin stiffness scaling, we 
have located T c for several couplings a. The results are 
graphed in Fig. [3] We compare our results with the ex- 
pression obtained by Liu?i& 



1 



1 



2 — cosk x cosk y + a(l — cosk z ) 



(5) 

We find that while this equation gives a reasonable es- 
timate for T c (a)/T c (l) for a close to unity, it begins to 
deviate substantially from the SSE results for small a. 



IV. CALCULATIONS OF THE SPECIFIC HEAT 

Having determined T c as a function of a, we now 
present the results for the specific heat calculations. The 
specific heat is defined as the temperature derivative of 
the energy, C v = (dE/dT)/N. As discussed in Appendix 
A, the SSE method allows us to obtain a direct estimate 
of the specific heat from the operator sequence in the 
simulation, so that any additional noise in the data due 
to numerical differentiation of the energy function can 
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FIG. 3: Ratio of the critical temperature for the anisotropic 
system to that for the isotropic system as a function of the 
anisotropy. The circles denote the results from Eq. 

be avoided (although the two approaches in practice give 
very similar results). The SSE estimator for the total 
specific heat (i.e., not normalized by the lattice size) is 

NC V = (n 2 ) - (n) 2 - (n), (6) 

where n is the power-series expansion order (the num- 
ber of bond operators in the SSE operator string), which 
fluctuates in the simulations. We will be interested in the 
contributions to C v from the spin-spin ordering across 
and within the layers close to T c . Decomposing the 
Hamiltonian into an in-plane term H p and an inter-layer 
term H z , the specific heat 

C v = (d(H p )/dT + d(H z )/dT)/N = CI + C*. (7) 

The SSE estimators for the two terms are given in terms 
of the numbers of bond operators in the expansion acting 
within a single layer (n p ) and between two layers (n z ): 

NCP = (n*) + {n p n z ) - {n p ) 2 - {n p ){n z } - (n p >, (8) 
NC* = (n 2 z ) + (n p n z }-{n z ) 2 -(n p ){n z )-{n z ).(9) 

These expressions suggest the possibility of a different 
decomposition of the specific heat. We will define 6 
as the part of the estimator © that contains only purely 
in-plane contributions: 

Cf ane = «n») - (n p ) 2 - (n p ))/N. (10) 

We refer to the remaining part of the total susceptibility 
as the 3D contribution, i.e., 

=c v - cp' anc = Q ntor + C™ oss , (11) 



FIG. 4: The specific heat over a wide range of temperature 
for several different anisotropies. The system size is 48 x 48 x 
12. The separation of the 3D ordering peak from the broad 
maximum arising out of the 2D physics is clearly visible for 
a < 2~ 3 . 

where the purely inter-plane contribution C^ ntcl and 
cross-term C£ ross are given by 

Q" tCr = {{<) - (n z ) 2 - (n z ))/N, (12) 
C c J oss = 2((n p n z )~(n p )(n z ))/N, (13) 

We will show that the cross-term, half of which appears 
both in Eq. JHJ and Eq. JJjJ, dominates in the 3D con- 
tribution (|llfl . The advantage of considering separately 
the different contributions to C v , either in the form of 
(J7J) or lfTT|) and, is that the full specific heat is dominated 
by the in-plane term and the other contributions can be 
difficult to discern due to statistical fluctuations. We will 
here focus in particular on the 3D contribution (fill) . 

The specific heat for the 3D Heisenberg model on 
highly anisotropic lattices (a < 1) will have two sepa- 
rate peaks, reflecting the 2D physics and the 3D ordering. 
The Mermin- Wagner theorem dictates that there can be 
no long-range order at T > in a strictly 2D system with 
a continuous symmetry. The correlation length then di- 
verges exponentially^ as T — > 0, and the specific heat has 
a broad maximum at T/J ss 0.7^ This broad maximum 
is the dominant feature of the specific heat curve also for 
small inter-planar couplings. On the other hand, for any 
a > there is a phase transition to an ordered state at 
T c > 0, as we have discussed in Sec. III. The signature of 
this phase transition in the specific heat should be a peak 
at T c . Since the transition belongs to the 3D Heisenberg 
universality class, there should a cusp- like singularity (in- 
stead of a divergent singularity) and the peak height is 
finite. 

SSE results for the specific heat over a wide temper- 
ature range are shown in Fig. 0] for a system of size 
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FIG. 5: The 3D contribution to the specific heat for several 
different anisotropies, a, and for three different system sizes. 
Results for the purely inter-plane term Ci" er are also shown 
for the three largest couplings (for the largest system size 
only). 



N = 48 x 48 x 12. The effects of finite system size on the 
position of the peak and the peak height will be discussed 
later. The separation of the 3D ordering peak from the 
broad maximum arising out of the 2D physics is clearly 
seen for a < 2~ 3 . It is also seen that the excess peak 
height over the 2D background decreases rapidly with 
decreasing a, becoming hard to discern for a < 2~ 5 . 

Since the specific heat curve is dominated by its 2D 
contribution when a <C 1, it is extremely difficult to 
study the nature of the 3D peak near T c . However, the 
3D contribution can be studied to a high degree 
of accuracy. Results for several couplings a and system 
sizes are shown in Fig. \5\ Several features are immedi- 
ately apparent. The 3D contribution peaks at the Neel 
temperature and rapidly decreases away from it. The 
peak position moves only slightly with increasing system 
size. The estimates of T c obtained from the position of 
the peaks are in close agreement with the more accurate 
estimates we obtained in Sec. Ill using the spin stiffness. 
In Fig. [5] we also show some results for the purely inter- 
plane contribution C™ ter to C 3D , which is seen to be 
small and decreasing relative to the full 3D contribution 
as a — > 0. This is expected, as the estimator <|12[) implic- 
itly contains a prefactor proportional to a 2 , whereas the 
cross-term (|13f) contains a linear a dependence. 

While the specific heat anomaly is most pronounced 
in the 3D contribution, it is also present in the purely 




FIG. 6: The specific heat and its in-plane contribution for 
a — 2 -4 . The anomalies at the transition temperature is 
clearly visible for both. The system size is 48 x 48 x 12. 
For comparison, the specific heat for the pure 2D Heisenberg 
model is also shown. 
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FIG. 7: The excess in C£ lano (T c ) over the specific heat of the 
pure 2D Heisenberg model at the 3D T c , normalized by the 
corresponding 3D contribution to the total specific heat. The 
system size is 48 x 48 x 12. 



in-plane term. This is shown in Fig. El where we have 
graphed the total specific heat and the purely in-plane 
contribution at a = 1/16, where the 3D ordering peak is 
well separated from the broad 2D maximum. We com- 
pare these results with the specific heat C™ for a 2D 
system (a — 0). As expected, the in-plane term for the 
3D system is dominated by a broad maximum and co- 
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incides closely with the specific heat of the 2D system 
away from T c . However, there is also a distinct peak at 
the 3D transition temperature. In order to quantify the 
relative sizes of the ordering peaks in C^ D and CP lanc , 
we next consider the excess at T c of the in-plane con- 
tribution over the specific heat of the pure 2D system 
model at the same temperature. Its ratio to the 3D con- 
tribution is graphed as a function of the coupling a in 
Fig. [3 As a — > 0, this ratio appears to converge to a 
value ~ 1, or, in other word, the ordering peak in the in- 
plane contribution becomes nearly equal to that of the 
3D contribution. 

The peak height C^ D (T C ) decreases rapidly with de- 
creasing a. To get a more quantitative estimate of the 
nature of its variation with a, we have extracted the ther- 
modynamic peak height for different a. The specific heat 
exponent, which governs the scaling of the peak to infi- 
nite size, is small (and negative)^ and the statistical 
errors of our data are relatively large for small a. The 
extrapolation is therefore affected by some uncertainty 
that is not easy to quantify precisely. Our results are 
shown in Fig. |H1 For small a, the peak height is nearly 
linear in a. This behavior can be roughly understood 
by the argument that the specific heat anomaly should 
scale as l/£ 2 , where £ is the correlation length of the 
2D system at the 3D transition temperature. Further- 
more, the 3D correlations become significant and lead to 
the 3D transition^ when £ 2 a rj 1. Thus the amplitude 
for the specific heat anomaly should vanish linearly with 
a. It would be interesting to compare the specific heat 
anomaly of various quasi-2D Hcisenberg systems against 
this result. 



V. CONCLUSIONS 

In this paper we have studied the 3D ordering tran- 
sition in a model of weakly coupled Heisenberg planes. 
Our results on the transition temperature and univer- 
sality class of the transition are in accord with general 
expectations. Our primary focus here was on the specific 
heat and in particular on the specific heat anomaly at 
the 3D ordering transition. We find that for small J± the 
amplitude for the specific heat anomaly is a nearly linear 
function of J± . It should be possible to compare this re- 
sult directly against experiments on various anisotropic 
materials. However, it is clear that for highly anisotropic 
systems (such as La2Cu04, where the anisotropy maybe 
as small as 10 -6 ) such anomalies will be very difficult to 
detect above the background. 
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FIG. 8: The peak height for the 3D ordering extrapolated to 
the thermodynamic limit as a function of the anisotropy. For 
small anisotropies, the peak height increases approximately 
linearly with the anisotropy. 
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APPENDIX A: THE SSE METHOD 

The SSE method has been discussed in several 
papers ^S°42i Here we present a brief outline of the 
method in order to discuss the estimator for the specific 
heat. For the present case, the SSE approach starts by 
casting the Hamiltonian in the form 

. 3N 

H = -~Y,[Hi,b-H 2>b ]+C, (Al) 

6=1 

where b denotes the bond connecting the nearest neigh- 
bor sites (i(b), j(b)) , C is an additive constant and the 
operators H\.b and H^b are defined as 

Hi,t = 2J{b)[\- S* i{h) S z m ], (A2) 
H 2 , b = J(b)[S+ b) S- {b) +S m S+ b) }. (A3) 

The coupling constant J(b) — J for bonds in the planes 
and J(b) = J± for inter-planar bonds. An exact and 
useful expression for an operator expectation value at 
inverse temperature /3 — J /T, 

(A) = -TriAe-P 6 }, Z = Tr{e-^}, (A4) 

is obtained by expanding the density matrix e~^ H in a 
Taylor series and writing the trace as sum over the diag- 
onal matrix elements in a basis {\a)} = {\Sf , . . . , Sfj)}. 
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The partition function can then be written as 

oo an n 

z = EEE ^rHlI^>>> ( A5 ) 

n=0 a S n p=l 

where S n denotes a sequence of index pairs defining the 
operator string J]p = i H ap . bp : 

S n = [Ol, 61] [03,62] ■ ■ ■ [On,K], (A7) 

where a G {1, 2, 3}, 6 G {1, ... , N}. We have separated 
the temperature dependence of the weight factor for con- 
venience. We can now write the expectation value of an 



operator as 

n=0 a S n 

Taking A = H, it can be showni2i24 that the energy is 
given by the average length of the operator sequences 

00 

' n=0 a S„ " 

A straightforward differentiation with respect to temper- 
ature gives the specific heat C v = |y in the form of 
Eq. p. 



1 L. L. Liu, and H. E. Stanley, Phys Rev. Lett 29, 927 (1972). 

2 L. J. de Jongh, and H. E. Stanley, Phys. Rev. Lett. 36, 
817 (1976). 

3 S. Chakravarty, Phys. Rev. Lett. 77, 4446 (1996). 

4 E. W. Carlson, D. Orgad, S. A. Kivelson, and V. J. Emery, 
Phys. Rev. B 62, 3422 (2000). 

5 M. Bocquet, F. H. L. Essler, A. M. Tsvelik, and A. O. 
Gogolin, Phys. Rev. B 64, 94425 (2001). 

6 I. Affleck, M. P. Gelfand, and R. R. P. Singh, J. Phys. 
A 27(22), 7313 (1994); I. Affleck, and B. I. Halperin, 
iMd.29(ll), 2627 (1996). 

7 S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys. 
Rev. Lett. 60, 1057 (1988); Phys. Rev. B 39, 2344 (1989). 

8 A. V. Chubukov, S. Sachdev, and J. Ye, Phys. Rev. B 49, 
11919 (1994). 

9 R. R. P. Singh, W. E. Pickett, D. W. Hone, and D. J. 
Scalapino, Comments on Modern Physics 2(1), Bl (2000). 

10 K. Sun, J. H. Cho, F. C. Chou, W. C. Lee, L. L. Miller, 
D. C. Johnston, Y. Hidaka, and T. Murakami, Phys. Rev. 
B 43, 239 (1991). 

11 G. S. Rushbrooke, G. A. Baker, and P. J. Wood, in Phase 
Transitions and Critical Phenomena, ed. C. Domb and M. 
S. Green, Academic Press, 1981. 

12 J. Oitmaa, C. J. Hamer, and Z. Weihong, Phys. Rev. B 
50, 3877 (1994). 

13 R. A. Sauerwein and M. J. De Oliveira, Mod. Phys. Lett. 
B 9, 619 (1995). 

14 A. W. Sandvik, Phys. Rev. Lett. 80, 5196 (1998). 

15 T. Oguchi, Phys. Rev. 133, A1098 (1964). 

16 S. H. Liu, J. Magn. Magn. Mater. 82, 294 (1989). 

17 V. Yu. Irkhin and A. A. Katanin, Phys. Rev. B 55, 12 318 
(1997); ibid.57, 379 (1998). 

18 I. G. Araiijo, J. R. de Sousa, and N. S. Branco, Physica A 
305, 585 (2002). 



19 A. W. Sandvik and J. Kurkijarvi, Phys. Rev. B 43, 5950 
(1991); A. W. Sandvik, Phys. Rev. B 56, 11678 (1997). 

20 A. W. Sandvik, Phys. Rev. B 59, R14157 (1999). 

21 O. F. Syljuasen and A. W. Sandvik, Phys. Rev. E 66, 
046701 (2002). 

22 E. Marinari, Lecture Notes in Physics, Vol. 501 Advances 
in computer simulation: lectures held at the Eotvos Sum- 
mer School in Budapest, Hungary, 16-20, July 1996, edited 
by J. Kertsz and I. Kondor (Springer, 1998). 

23 K. Hukushima, H. Takayama, K. Nemoto, Int. J. Mod. 
Phys. C 7, 337 (1996); K. Hukushima, K. Nemoto, J. Phys. 
Soc. Jpn. 65, 1604 (1996) 

24 P. Sengupta, A. W. Sandvik, and D. K. Campbell, Phys. 
Rev. B 65, 155113 (2002). 

25 W. Kohn, Phys. Rev. 133, A171 (1964). 

26 P. Kopietz, Phys. Rev. B 57, 7829 (1998). 

27 E. L. Pollock and D. M. Ceperley, Phys. Rev. B 36, 8343 
(1987). 

28 K. Harada and N. Kawashima Phys. Rev. B 55, R11949 
(1998). 

29 A. Cuc coli, T. Roscilde, V . Tognetti, R. Vaia, and P. Ver- 
rucchi, cond-mat/0209316 

30 A. W. Sandvik, Phys. Rev. Lett. 83, 3069 (1999). 

31 For a review of finite-size scaling, see: M. N. Barber in 
Phase Transitions and Critical Phenomena, Vol. 8, ed. 
Domb and Lebowitz (Academic Press, 1983). 

32 M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi, 
and E. Vicari, Phys. Rev. B 65, 144520 (2002). 

33 G. Gomez-Santos, J. D. Joannopoulos, and J. W. Negele, 
Phys. Rev. B 39, 4435 (1989); J.-K. Kim and M. Troyer, 
Phys. Rev. Lett. 80, 2705 (1998). 

34 D. C. Handscomb, Proc. Cambridge Philos. Soc. 58, 594 
(1962); 60, 115 (1964). 



